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Abstract 

We present a method for solving the reaction-diffusion equation with general potential in free space. It is based on 
the approximation of the Feynman-Kac formula by a sequence of convolutions on sequentially diminishing grids. 
For computation of the convolutions we propose a fast algorithm based on the low-rank approximation of the Hankel 
matrices. The algorithm has complexity of 0(nrM log M + nP'M) flops and requires 0{Mr) floating-point numbers 
in memory, where n is the dimension of the integral, r <s: n, and M is the mesh size in one dimension. The presented 
technique can be generalized to the higher-order diffusion processes. 
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1. Introduction 


Path integrals dciia play a dominant role in description of a wide range of problems in physics and mathematics. 
They are a universal and powerful tool for condensed matter and high-energy physics, theory of stochastic processes 
and parabolic differential equations, financial mathematics, quantum chemistry and many others. Different theoretical 
and numerical approaches have been developed for their computation, such as the perturbation theory ||4l, the station¬ 
ary phase approximation 13161, the functional renormalization group QEl, various Monte Carlo m and sparse grids 
methods IfTOlfTTI . The interested reader can And particular details in the original reviews and books IfT^[131 fT4l . 

In this paper we focus on the one-dimensional reaction-diffusion equation with initial distribution f(x) : K ^ 
and a constant diffusion coefficient cr 


— m(x, f) = cr — rM(x, f) — y(x, t)u(x, t), 
ot dx-^ 

u(x, 0) = fix) 


te[0,T], xeM. 


( 1 ) 


This equation may be treated in terms of a Brownian particle motion lITSl [161 [T71l . where the solution m(x, f) : K x 
[0, T] is the density distribution of the particles. The potential (or the dissipation rate) V(x, t) is bounded from 

below. We do not consider the drift term p-^u(x,t) because it can be easily excluded by a substitution u{x,t) —> 
m(x, ITSl . 

The solution of Q can be expressed by the Feynman-Kac formula IfTSi [T9l l20l 


u(x,T)= f 

Jcix,0;ri 


( 2 ) 


where the integration is done over the set C{x, 0; T) of all continuous paths ^(T) : [0, T] ^ M from the Banach space 
H([0, T],®) starting at ^(0) = x and stopping at arbitrary endpoints at time T. is the Wiener measure, and ^(f) 


Email address: m. litsarevOskoltech.ru (Mikhail S. Litsarev) 


Preprint submitted to Journal of Computational Physics 


November 9, 2015 






is the Wiener process One of the advantages of the formulation (|^ is that it can be directly applied for the 

unbounded domain without any additional (artificial) boundary conditions. 

Path integral (|^ corresponding to the Wiener process is typically approximated by a finite multidimensional 
integral with the Gaussian measure (details are given in Section [ZT] ). The main drawback is that this integral is a high¬ 
dimensional one and its computation requires a special treatment. Several approaches have been developed to compute 
the multidimensional integrals efficiently. The sparse grid method l24l has been applied to the computation of 
path integrals in ll25l . but only for dimensions up to ~ 100, which is not enough in some applications. The main 
disadvantage of the Monte Carlo simulation is that it does not allow to achieve a high accuracy EglETl for some 
cases (highly oscillatory functions, functions of sum of all arguments). 

The multidimensional integrand can be represented numerically as a multidimensional array (a tensor), which 
contains values of a multivariate function on a fine uniform grid. For the last decades several approaches have been 
developed to efficiently work with tensors. They are based on the idea of separation of variables ll^l^l30ll3Tll firstly 
introduced in It allows to present a tensor in the low-rank or low-parametric format llT4l[T5lf3^ . where the 

number of parameters used for the approximation is almost linear (with respect to dimensionality). To construct such 
decompositions numerically the very efficient algorithms have been developed recently: two-dimensional incomplete 
cross approximatiot^ for the skeleton decomposition, three-dimensional cross approximation IfJTl for the Tucker 
format ll38lf39ll40ll4lT in 3D, tt-cross i42l approximation for the tensor train decomposition MM, which can be also 
considered as a particular case of the hierarchical Tucker format ll45ll46ll47l for higher dimensional case. For certain 
classes of functions commonly used in the computational physics (multiparticle Schrodinger operator ll48l 149115011511 
|^|53, functions of a discrete elliptic operator Il5^l55ll5^l57ll58ll59l . Yukawa, Helmholtz and Newton potentials ll60l 
|6l]|62]|63l, etc.) there exist low-parametric representations in separated formats and explicit algorithms Il64ll65l to 
obtain and effectively work with them (especially quantized tensor train (QTT) format 11^1^7116811^1701171117211771 1. 
In many cases it is very effective to compute the multidimensional integrals m using separated representations GS, 
particularly for multidimensional convolutions ITT^ ITTI l78l 1791 and highly oscillatory functions llSOl . 

Our approach presented here is based on the low-rank approximation of matrices used in an essentially different 
manner. We formulate the Feyman-Kac formula as an iterative sequence of convolutions defined on grids of dimin¬ 
ishing sizes. This is done in Section 3.2 To reduce the complexity of this computation, in Section [33] we find a 


low-rank basis set by applying the cross approximation (see Appendix A i to a matrix constructed from the values of 
a one-dimensional function on a very large grid. That gives reduce of computational time and memory requirements, 
resulting in fast and efficient algorithm presented in Section [34| The numerical examples are considered in Section]^ 
The most interesting part is that we are able to treat non-periodic potentials without any artificial boundary conditions 
(Section |4.3| l. 


2. Problem statement 

2.7. Time discretization 

Equation Q corresponds to the Wiener process. A standard way to discretize the path integral is to break the time 
range [0, T] into n intervals by points 


Ti^-k-6t, 0 < k < n, n: T„- T. 

The average path of a Brownian particle firf) after k steps is defined as 

^f{Tk)^X + ^l+f2 + ...+^U, 

where every random step 1 < i < k,is independently taken from a normal distribution 7V(0, 2cr6t) with zero mean 
and variance equal to 2cr6t. By definition, - x. 


'Because the low-rank representation of large matrices based on the adaptive cross approximation is directly related to the manuscript we 
summarize the basics of the method in|Appendix A| 
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Application of a suitable quadrature rule on the uniform grid (i.e., trapezoidal or Simpson rules) with the weights 
to the time integration in (|^ gives 


A(T)^ f V{^{T),T-T)dT^^WiVf6U = 

"^0 i =0 


and transforms the exponential factor to the approximate expression 


e 


-Air) 


O' 

1=0 


ViV'"^St 


The Wiener measure, in turn, transforms to the ordinary n-dimensional measure 


( 3 ) 


£)?’ = I - 




4crdf ’ 


and the problem reduces to an n-dimensional integral over the Cartesian coordinate space. Thus, we can approximate 
the exact solution (|^ by f) 

u(x, T) - lim T) 

written in the following compact form 


X oo ^ 

“ i=0 


(4) 


The integration sign here denotes an n-dimensional integration over the particle steps and is defined in Q. 
The convergence criterion in terms of n for the sequence 0 is discussed and proven in 1171 . p. 33. The limit of Q 
exists if it is a Cauchy sequence. 

Our goal is to compute the integral Q numerically in an efficient way. 


3. Computational technique 

3.1. Notations 

In this paper vectors (written in columns) are denoted by boldface lowercase letters, e.g., a, matrices are denoted 
by boldface capital letters, e.g., A. The i-th element of a vector a is denoted by a,, the element (/, j) of a matrix A 
is denoted by A,y. A set of vectors a^, niQ < m < mi is denoted by {am)mLm„, and the /th element of a vector a,„ is 
denoted by a„,i- 


Definition 1. Let a e and b e K'” be vectors and k > m. We say that a vector c e is a convolution of two 

ordered vectors a and b and write 

c = a o b, 

if c has the following components 


m-l 

Cl — ^ ai+jbj, ai — 0, V/: {/1 i < 0 V / >= k}. 
j=o 

The computation of the convolution can be naturally carried out as a multiplication by the Hankel matrix. 

Definition 2. We say that the Hankel matrix A 6 is generated by row a^ 6 and column b 6 and denote 
this by 

A-[a^b]H, 
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' ao 

ai 

02 ■ 

ak-2 

Ok-f 

Ql 

02 

03 

Ok-l 

bo 

02 

03 

04 ■ 

bo 

bi 

ak~2 

Ok-l 

bo ■ 

bk-4 

bk-i 

,ak-\ 

bo 

b\ ■ 

bk-3 

bk-^2, 


o’- - {aQ,a\, 02, ■■■,ak-2,ak-i), ^ - {bo,b\,b2, ■ ■ ■ ,bk~2)^ ■ (5) 

This compact notation will be used to compute convolutions (when they are written as a Hankel matrix-vector prod¬ 
ucts). As it can be directly verified, Va e K 

a ■ A. - [a ■ a^, a ■ b]//- (6) 


Definition 3. For two vectors a and b from Q for the case a, - bi,'ii : Q < i < k - I, we will also write 



This notation will be used when vector b is a subvector of a. 


(7) 


3.2. Multidimensional integration via the sequence of one-dimensional convolutions 

Multidimensional integral Q can be represented in terms of an iterative sequence of one-dimensional convolu¬ 
tions. Indeed, for a one-dimensional function such that 





^+1 




X e. 


k — n,n — 


(8) 


with 

and the initial condition 
the solution Q reads 


= fix), 


(9) 

( 10 ) 

( 11 ) 


The iteration starts from k = n and goes down to k - 1. Since the function <l)^"\x) is bounded and the convolution (j^ 
contains the exponentially decaying Gaussian, the integral has finite lower and upper bounds. Consider 


Ff\x) ^ Ffix) ^ ^ + f) (12) 

We suppose that the product <l>^^\(.r H- f)e~'^^ rapidly decays, so that for a^ large enough, we can approximate the 
integral F^"\x) in ([^ by F^^\x) and assume that this approximation has an error e in some norm 

\\Ff\x)-Ff\x)\\<e. 


This approximation has an important drawback; as soon as Ff\x) has to be computed on the semi-open interval 
[-Ojc, Ox), the domain of fI”\x) should be taken larger, i.e. [-na^, naf) for n steps, because of the convolution structure 
of the integral ([T^. Indeed, if we suppose, that the function F^^\x) is computed on the uniform mesh 


X® = —kox -H ihx. 


M = 2Nx, 


0 < i < kM, 
4 


hx — OxINx, 


(13) 




^0 ^10 ^20 ^0 ^10 fto 



Figure 1: A correspondence of meshes for two nearest iterations from equation for k = 1. Blue filled circles sepai'ate the 

ranges coiTesponding to dilferent steps m, 1 < m < k in time [-max, max). Ticks on the axes label the mesh points. Violet curved lines show 
correspondence CD between the two meshes for nearest iterations. 


and the integration mesh is chosen to be nested in ( [T3] ) with the same step 

^ -a^ + jh„ 0<j<M, (14) 

then the function (x) is defined on the mesh 

xf^^'>^-(k+l)a^ + ih^, 0<i<{k+\)M, (15) 


and 




(16) 


The last equality follows from definitions ( [T3] l and ( [T4l l. This is illustrated in Figure[2 

The integral ( [T^ can be calculated for every fixed xf^ of the mesh ( [T3| ) as the quadrature sum with the weights 




M-\ 


7(n) 


[xf) ^ Yj ^i\ Pi'l. 

1=0 




(17) 




(18) 


The complexity of the computation of A"' O(kN^) flops. It can be reduced to 0{kNxlogNx) by ap¬ 

plying the Fast Fourier transform (FFT) for convolution ( fTT] ). Full computation of Pf'' (jc-*') costs 0{rPNx\ogNx) 
operations and 0{nNx) floating-point numbers. This complexity becomes prohibitive for large n (i.e., for small 
time steps), but can be reduced. Below we present a fast approximate method for the calculation of F*"* in 
0{nrNx\ogNx + nr^Nx) flops and 0{rNx) memory cost with r n, by applying low-rank decompositions. 


3.3. Low-rank basis set for the convolution array. 

In this section we provide a theoretical justification for our approach. Consider a sequence of matrices 6 
■^kMxM corresponding to the iterative process © and constructed in the following way 

A® = (x®.), 0<j<M, 0<i<kM, (19) 

where k is the iteration number. 
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Figure 2: Transition between two neighbour iterations is illustrated. The left-hand-side matrix A is multiplied by vector p in a resulting vector s 
according to j20^ . Explicit structure of matrix blocks H„, in and vector blocks in 0 is shown. Then entries of vector s are multiplied by 
corresponding factor to produce the next iteration step From a new vector a' there formed a new matrix A' according to 0. 

The last point from the previous iteration is not needed and is thrown out. Then the steps repeat for the next iteration. 


Let us now consider iteration ( [T7| i at the step k - and for simplicity omit the index A:o + 1 in the matrix and mesh 
notations in- Let us also denote the sum for Xi taken from the grid ( [T3) by i,- = (x,) and set pj = Pjp(A, ^j). 

Then 

M-l 

Si = ^ Aij Pj « s = Ap. (20) 

2=0 

The equality ( |20l l establishes the recurrence relation between iterations at the step k (the right-hand side) and the 
step k - \ (the left-hand side) according to ( [T7] i and ( [TSl l, see Figure]^ 

The matrix A is a Hankel matrix, as follows from definition ( [T9] l, and consists of k square blocks H^, 0 < m < k, 
such that 

A^ = (Ho, Hi... Hu). (21) 

Here, every block H,„ is a Hankel matrix as well generated by the upper row 1,^ and the right column r„,+i correspond¬ 
ingly; 

H,„ = |^1„,, r„,+i , 

(the notation [a, b]// is introduced in Section [TT] Dehnition]^, where 


- («io> Oio+1> • • ■ ^ aio+M-2, O/o+M-l); 
Lk - («io> Oio+l, . . . , ai^+M-l), 


io — m ■ M, 


0 < m < k. 


( 22 ) 


by the definition, see Figure]^ It can also be represented as a sum of two anti-triangulaij^ Hankel matrices 

H,„ = L,„ -I- R„,, Lm = r(m+l) M]^ . (23) 

where the upper-left L,„ has nonzero anti-diagonal and the bottom-right R„, has zero anti-diagonal according to ( |22| l. 
Equation ( |20l l may be rewritten in the block form (see again Figure]^ 


h,„ — H;jjP, 




( 24 ) 


^By anti-triangular matrix we call a matrix which is triangular with respect to the anti-diagonal of the matrix. 
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Here every block H„ is multiplied by the same vector p. The number of matrix-vector multiplications can be reduced, 
if the dimension d of the linear span “N = is less than k. Before estimation of the dimension we formulate 

some auxiliary lemmas (proven in the [Appendix C| l. 

Lemma 1. Let be a basis set of span r\ < k, and let matrix U,- = oj^. Then {U/lJlg* is a basis set 

of span {L„)^;i'g/rom 

Lemma 2. Let {w/ljig* be a basis set of span ^2 < k, and let matrix W; = |^0^, w,j^. Then is a basis 

set of span from ( |2^ . 

Lemma 3. Let be a basis set of span — ('^J y “/,(m- 1 )) according to Q- Then {w/)^^g* is 

a basis set of span {fm)* 

Let us define a basis set * as follows 


Q,= 



0 <i < r 
r <i <2r 


(25) 


An obvious corollary of the previous Lemma is the following Theorem. 


Theorem 1. The dimension of the linear span of matrices is equal to 2r. Moreover, it is contained in the 

linear span of the matrices {Q/)^i^o * defined in 


Proof The matrix can be written as a sum ( |2^ , = L,„ + R,„. According to Lemma^ the set {U,)^^q is a basis 

set of the span {Lml^^Lg. By Lemma 5 {WillTg is a basis set of the span and by Lemma 2 {WIJTa 

_ t. 1 ^ _- ^ 1 I—I - I—I 


setof{R)-Jg. 
is 2r. 


The subspaces {UijJLo and {W, 


'.)&1 


is a basis 

contain only zero matrix in common, so the dimension of the basis 

□ 


Lemma 4. Let {u/)J'^q be a basis set of span {lm)m=o- Then for basis matrices {Qil^^Tg * defined in P5| l the computation 
of the matrix-by-vector products 

k, = U,p, t, = W,p, (26) 

costs 0{M log M) flops for a fixed 0 < i < r. 


Proof Consider a Hankel matrix 



A product G,p is a result of the convolution u, o p, which can be done by the FFT ISTl |82l procedure in 0{M log M) 
flops for a fixed 0 < i < r. The vector p = (pm-i, ■ ■ ■ ,,Pi, PoV is taken in the reverse order. □ 


Once the basis {Q()f,,o* found, the complexity of the multiplication Ap in ( |20l l can be estimated 

as follows. 

Theorem 2. Let the set {Q, )flo^ defined in ( |25| ) be a basis set of the linear span '7L generated by the set of Hankel 
matrices H,„ defined in •ID- Then the computation of any elements s, of the vector s ( |20l i costs 0{rM log M r^M) 
flops for Ks — O(rM). 

Proof. Indeed, by the assumption H,„ = 2^=0^ Cm,Q, for each m, 0 < m < k. The complexity of the product Q,p, 
0 < i <2r for a fixed i is 0(M log M) flops by Lemma|^ The computation of such products for all i takes OirM log M) 
flops. 

The vector h„„ which is a subvector of s, is represented via few matrix-by-vector products ( |26] l as follows 

r r 

kff) — H,nP — L^p R;;;P — ^ tTm/U/P ^ ] O^ffi/k/ + (27) 

1=0 1=0 


The computation of its /th component h^i takes 0{r) flops for any m. Computation of O(rM) components sj of the 
vector s, which are also the components of the particular vector h^ (for m = [ j^J), costs, in turn, O(r^M) flops. Finally, 
C)(rM log M-H r^M). □ 
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Remark 1. Each component of the resulting vector can be computed by the formula 




1=0 


LM 


Ij - i mod M. 


(28) 


Here kn^ is the Ij-th component of the vector k, and tn. is the l^-th component of the vector t, . 
Remark 2. It follows from Lemma[^that cr,+ij = /?,) in ( |27l i. 


3.4. Final algorithm 

To compute F^"\xi), which defines the final solution ( [TT] i on the mesh ( pj) , one needs to carry out iterations ( [T7] i 
starting from k - n down to 1. At each iteration step k we construct a function fk which approximates the 
entries s] in equation ( |28] l as follows. Suppose, that the function fk+\ has been already constructed at the 

previous step A: + 1. Then, to compute fk at the current iteration k, we consider the matrix with the entries 






-wicV(yij,T„-t)6t 


yij = 


_ j^+i) 




(29) 


and apply the cross approximation ( |A.2[ ) to this matrix. The columns of this matrix are vectors element-wise 
multiplied by the corresponding exponential factor with the potential ( |29l l, see Figure The algorithm of the cross 
approximation requires only O(rM) entries, which are being chosen adaptively. They are calculated by the function 
fk+\ (yij^ on-the-fly for the particular points yij. Thus, 


^(^+1) ^ Bv^, B 6 V e r^M, 


(30) 


where B and V are matrices of the rank r saved in memory. By construction, the m-th column of matrix is the 

vector Im from ( |2^ and the /-th column of matrix B is the basis vector u, from Lemma [I] Hence, is the matrix of 
coefficients of the decomposition ( |C.l[ l. Once the cross approximation ( [30l l is obtained, the memory allocated for all 
data structures related to fk+\ (x(*^*'jcan be overwritten at the next iteration. 

Computation of the circulant matrix-vector products ( |26| l is done according to Lemma |^by the convolution g,- = 
b, o p, where b, is a column of the matrix B. The vectors g, = (t,, k,)^ are also saved in the memory. Then fk is 
calculated by equation ( |28l ), and the algorithm proceeds to the next iteration. 

At some iteration step k the rank of the decomposition ( [30l l will reach the number of columns and from this 
iteration it will be more efficient to carry out the convolution (20i without low-rank approximation. Complexity of 
one iteration of the presented algorithm is estimated in Theorern|2 Finally, for all n steps it is 0{nrM log M + nr^M) 
flops, r ^ n. The standard LET based algorithm applied to the whole array without any low-rank compression at each 
step gives complexity for all n steps equal to 0{n^M log M) flops. We illustrate this theoretical estimations by the 
example from Section [42| in Figure]^ 

Basically, the asymptotic complexity, proven in Theorem]^ is practically useful for r <s: n. This is the main as¬ 
sumption for the matrix from ( [29| ). Existence of such an approximation (and the properties of the initial problem) is in 
general still an open question. Some particular cases were studied in If83l . It was shown, that the cross approximation 
converges for matrices having singular vectors satisfying the coherence property. Some estimations can be found 
in Il84l[85l[8^ also. There is a theoretical idea how to identify the existence of the low-rank structure of a given matrix 
generated by a one-dimensional target function a priory (see lISTll and Appendix B for details). 


4. Numerical experiments and discussions 


4.7. Flarmonic Oscillator 

As a first example, let us consider a model system, which can be solved analytically, with the initial condition 
fho(x) and the dissipation rate Vhoix, t) defined as 


fho(x) = p(J3, x) = 



^ho{x^ t) 


t+l' 


(31) 


^but do not compute all its elements 


8 












-IQx 

-6ax ■ 


-Aux 

-3ax 

-2ax 


0 


Gx 

IGx 

3gx 

AGx 


6gx 7a 

t 

t 


.t. 

" " f ' 

11 Ml 111 

.t. 


t 

1 Ml 111 

. t "" 

t 

Mill 

.t '" 

t 

1 

^(7) 

1 

J.1) 


1 

v(7) 

T 

AD 


1 

AD 


1 

x(D 


1 

x(D 

1 

J.D 


1 

<j) 

1 

AD 

^0 



-^30 

''■45 


■^60 


4,75 


XgQ 

'^'105 


''•120 

^135 
















■'o 

i 

■^15 

i 

11 1 1 #! 1 11111 

1 I#! 11 1 

■'30 

i 

1111 I#! 1111 

■*45 

i 

1 1 1 11 1 1 1 1 1 

MINI 

i 

1 1 1 1 1 1 1 1 1 

mm- 

3,5 

i 

mm 

mil 

■*90 

i 

1 1 1 1 1 1 1 1 

■'105 

i 

1 1 1 1 1 1 1 1 1 1 1 

1 1 1 1 1 

■'120 

i 

1111 1 111 

■*135 

i 

1 1 1 1 1 1 1 1 1 1 1 1 1 



Qo 

aio 

O40 

060 

080 

O|00 

O12O 


Oo.o 

Ol,0 

02,0 


vo,o 

Vo,l 

Vo,2 

■ ■ Vo,17 4 'o,18 Vo,19 

a\ 

021 

O41 

061 

081 

OlOl 

O12I 


00,1 

Ol.l 

02,1 


Vi,o 

Vl.l 

Vl,2 

• ■ Vi,i 7 Vi_i 8 Vl,l 9 

0.1 

022 

042 

062 

082 

O|02 

O122 


Oo,2 

Ol,2 

02,2 


V2,0 

V2,l 

V2,2 

■ ■ V2,17 4'2,18 V2,19 















v^ 

an 

On 

O57 

O77 

O97 

0117 

O137 


Oo,17 

Ol,17 

02,17 



go 

= U() 0 p 

ai8 

038 

058 

078 

O98 

O1I8 

0138 


O0,18 

Ol,18 

02,18 



gl 

= Ui 0 p 

ai 9 

039 

O59 

O79 

Ggg 

0|19 

Ol39 


Oo,19 

Ol,19 

02,19 

B 


g2 

= U2 0 p 


Figure 3: Construction of matrix from a one-dimensional convolution l |20) according to algorithm in Section]^] On a spatial homogeneous 
mesh O the corresponding entries of vector s are calculated. By definition, vector s is composed from vectors h,„ ID- Each column of 
the matrix is composed of h,„ multiplied by a corresponding factor Then this matrix is decomposed by a cross approximation 

^{k) _ gyr por the approximation there needed only some elements of matrix which are chosen adaptively and computed on-the-fly. 
Then convolutions g, = u, o p ai'e calculated via fast Fourier transform and saved in the memory. Pai'ticular values of for the next iteration step 
k - \ can be computed by formula then. 


CPU Time 



Dimension n 


Figure 4: A numerical illustration of theoretical estimations for the example from SectionjA^ For a standard FFT based algorithm applied to the 
whole ari'ay the Oin^MlogM) flops complexity is labeled by square points. The low-rank complexity 0(nrM log M + nr^M) flops is labeled by 
circles. The time is scaled in minutes, n is the number of dimensions (iteration steps), M = 8000, r = 10. 
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According to equation ([TT]i the exact solution t) for the particular case (|3T]) has the following form (see Ap¬ 


pendix D for derivation) 


Comparison of the numerical low-rank solution with the exact one (|32l) gives the relative error 


(32) 


u-u 


(33) 


which in the order of magnitude is equal to the machine precision, where u is an approximate solution on the final 
mesh and u is the exact one on the same mesh. For our example 


Hi 




{x,) e 


-woVi,o(xi,T)Sr 




-«oVho{Xi,T)6t 


Here cr - 0.25, T = 10, n = 100, and the mesh is a uniform one on [-2,2] with M - 2Nx — 8000 points. It is 
interesting that the scheme is exact for this case. 


4.2. Cauchy Distribution 

The second example is taken from ll25l and is interesting because it can be solved analytically as well. For Vdx, f) 
and initial condition fdx) such that 


Vdx, f) = - 


1 

t + 1 


+ 2cr 


3x2- j 
(x2 + 1)2 ’ 


fdx) = 


1 1 
TT x2 H- 1 ’ 


(34) 


the exact solution is 


. ^ f 

udx, t) — — 


n x^ + \ 


In TabIe[T]we present numerical results demonstrating the numerical order of scheme by the Runge formula 


P = log2 


||U„ - U„/2|| 


with respect to 6t and the timings for the whole computation. Here u„ is the computed solution at the final step in 
time. 

Using our approach, it becomes possible to calculate u^’'\x, t) for large values of final time T due to th e low- 
rank approximation of matrices composed from the columns of the integrand values (see Section 3.4 1 . That 


significantly reduces the computational cost. For an example, for the last row of Table [T] iterations start from the 
calculation of the convolution on the range [-16386,16386) with 32 772 000 mesh points. This is reduced to the 
calculation of 10 (the rank) convolutions of two arrays with 8000 elements. 

As it can be seen from our results, the scheme has the second order in time. It can be improved to higher orders 
by Richardson extrapolation on fine meshes Il88] [89l . Another way is to use other path integral formulations with 
high-order propagators IIOOllTO . 


4.3. Nonperiodic potential with impurity 

The dissipation rate V{x, t) causes the creation and annihilation of diffusing particles, as it follows from the main 
equation Q- Without the Laplacian, which is responsible for the free diffusion, we have 

—u(x, t) - -V(x, t)u(x, t). 
at 

It can be seen, that the density of particles increases over time for V{x,t) < 0 and decreases for V(x,t) > 0 corre¬ 
spondingly. The case V{x, f) < 0 may lead to an instability in the solution, because the integral 


X CO 

fix + 

oo 


(35) 


10 










Table 1: Convergence rate for system H). Accuracy of the cross approximation £ = 10 Direct convolutions start from n = 20, cr = 0.5, range 
of final spatial domain is [-2,2), Nx = 4000. Dimension of the integral 0 is labeled by n, 6t is a time step, T is a final time for solution u{x, T), 
6 is an error estimated by the Richardson extrapolation, and p is the order of the scheme for 6t. Ranks of the matrix from (H are presented 
in column labeled by r. The CPU time for computation of the integral in all points of the mesh is reported in the last column. 


T 

n 

6t 

P 


e 

r CPU Time (min.) 

1.0 

32 

3.1 

10-2 

- 

2.8 


10 

0.1 


64 

1.6 

10-2 

- 

7.0 

10-2 

10 

0.2 


128 

7.8 

10-2 

1.997 

1.8 

10-2 

10 

0.4 


256 

3.9 

10-2 

1.999 

4.4 

10-® 

10 

0.9 


512 

2.0 

10-2 

2.0 

1.1 

10-® 

10 

1.8 


1024 

9.8 

10-4 

2.0 

2.8 

10-2 

10 

3.8 

20.0 

32 

6.3 

10-^ 

- 

4.1 

10-2 

10 

0.1 


64 

3.1 

10-2 

- 

1.6 

10-2 

10 

0.2 


128 

1.6 

10-2 

1.10 

4.8 

10-2 

10 

0.4 


256 

7.8 

10-2 

1.68 

1.2 

10-2 

10 

0.9 


512 

3.9 

10-2 

1.93 

3.1 

10-2 

10 

1.9 


1024 

2.0 

10-2 

1.98 

7.9 

10-4 

10 

4.0 


2048 

9.8 

10-2 

1.995 

2.0 

10-4 

10 

8.0 


4096 

4.9 

10-2 

1.999 

4.9 

10-2 

10 

16.8 


8192 

2.4 

10-2 

2.0 

1.2 

10-2 

10 

37.5 



Figure 5: Potential V{x) and initial distribution f{x) for periodic system with impurity ^36^ . Potential oscillates on a free space. Functions V(x) 
and f{x) are relatively shifted to break the symmetry. 
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u(x,t) 



Figure 6: Convergence of solution u{x, t) for nonperiodic potential with impurity j36) for different n. This results correspond to the data presented 
in Table The number of spacial mesh points M = 2Nx = 8000 in the final range [-2,2). The dissipation rate i |36) leads to a decrease in the norm 
of the distribution density. As seen in the picture, the solution is far from the correct one for the dimensions n = 64,128,256. 


Singular values 



0 200 400 600 800 1000 


Iteration step 


Figure 7: The first few singular values (s.v.) of the matrix for system 0 at each iteration step. The first s.v. cri is presented in the absolute 
value. The other ones are given in the relative values as (TiIcti . The values below the cross accuracy s = 10“^^ are thrown out. As it can be seen, 
approximate SVD-rank is similar to the cross rank (in the sense of criterion ^A.3)). 


12 












Table 2: Convergence rate for system j36) . Accuracy of the cross approximation e = 10 Direct convolutions start from n = 20, cr = 0.25, final 
domain is [-2,2), Nx = 8000. Dimension of the integral a is labeled by n, 6 t is the time step, T = 20 is the final time. The order of the scheme 
p 2 for 6 t and the relative error ei are estimated from the original data computed by the algorithm from Section [T^ The next values p 4 and €4 
are estimated by the Richardson extrapolation. As it can be seen, the scheme has the fourth order in time after the extrapolation. The ranks of the 
matrix fromi |30^ are given in the column labeled by r. The CPU time for computation of the integral in all points of the mesh is reported 
in the last column._ 


n 


6t 

Pi 


^2 

P4 


^4 

r 

CPU Time (min.) 

64 

3.1 

■ 10 -^ 

- 


- 

- 


- 

9 

0.2 

128 

1.6 

■ io-‘ 

- 

8.3 

■ 10-2 

- 


- 

9 

0.3 

256 

7.8 

■ 10-2 

1.47 

3.3 

■ 10-2 

- 

2.8 

■ 10-2 

9 

0.8 

512 

3.9 

■ 10-2 

1.62 

1.1 

■ 10-2 

2.00 

7.0 

■ 10-4 

9 

1.7 

1024 

2.0 

■ 10-2 

1.84 

3.1 

■ 10-2 

3.04 

8.6 

■ 10-2 

9 

3.6 

2048 

9.8 

■ 10-2 

1.95 

8.1 

■ 10-4 

3.66 

6.8 

■ 10 -'’ 

9 

7.0 

4096 

4.9 

■ 10-2 

1.988 

2.0 

■ 10-4 

3.85 

4.7 

■ 10-2 

9 

14.7 

8192 

2.4 

■ 10-2 

1.997 

5.1 

■ 10-2 

3.98 

3.0 

■ 10 -* 

9 

33.0 


may diverge (see Eq. 0 )- Therefore, when choosing V(x, t) < 0, one should make sure that the integral in ( [T5] l 
converges. 

Consider the following problem (see Figure]^ 

Vi(x) - a + sin^ (nl- + li)-- 3 -, fi(x) = a - 0.5, (3 = 0.5. (36) 

^ // + V;r 

It can be interpreted as a nonperiodic system with an impurity. The term V(x) does not decay in the spatial domain 
and it is not periodic. Therefore the reduction of this problem to a bounded domain is not a trivial task and would 
require sophisticated artificial boundary conditions. 

In Tablej^we present results of numerical calculations, which show the order of the numerical scheme. In Figure]^ 
we also present the computed solutions for different values of n. Even in this case, the solution converges with the 
order p - 2. We also used the Richardson extrapolation of u(x, T) for different n to get higher order schemes in time. 

4.4. Monte Carlo experiments 

In this section we present results of Monte Carlo simulation. To estimate the solution in a fixed point xq the 
following formula is used 

k=l i=0 

^(k){i) - ^(k)\ + ... + ^{k)i, 

where each component of the vector ^(k) - {^(k)i, ■ ■ ■ ,^(k)nY is independently taken from the normal distribution 
N(0,2cr6t) at each trial step k: \ <k < K, where K being the number of trials. 

Results for the exactly solvable model are presented in Table[^ We compare accuracy and timings for Monte 
Carlo and low-rank calculations. It should be emphasized that in the Monte Carlo approach only one point of u{xq, T) 
is calculated for a fixed xq in one simulation, while our approach allows to compute the whole array u{xi, T) on the 
whole mesh simultaneously. This numerical experiments have been done on a single CPU core without parallelization 
of the Monte Carlo algorithm just to estimate the speedup of the low-rank computation. More advanced realization 
such as quasi Monte Carlo methods can be used. As it can be seen, the low-rank algorithm presented in Section|30is 
much faster. 


5. Conclusion and future work 

The presented results show that the proposed method is an efficient approach for solving diffusion equations in a 
free space without artificial boundary conditions (ABC). Instead of standard solvers based on the ABC designed for 
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Table 3: Timings for system (34). Accuracy of the cross approximation e = 10 Direct convolutions start from n = 30, cr = 0.5, range of final 
spatial domain is [-2,2), Nx = 4000. Dimension of the integral is labeled by n, 6t is a time step, T is a final time for solution u{xo, T) computed 
in a fixed point jcq- Here xq = 0, T = 1. The relative en'or 6 = \u{xo, T) - w(xo, T)\/\u(xo, r)| is computed in one point xq. Time for one point 
calculation is presented for Monte Carlo approach ilz) and is estimated for the whole mesh array consisting of M = 2Nx = 8000 points (the last 
column). For the low-rank computation the total timings are presented as well. Monte Caido simulation has been done with K = 10^ samples. The 
low-rank results are labeled by LR, while the Monte Caido results are labeled by MC. 


n 

6 t 

m(xo, T) 


e 

CPU Time (\ point) 

CPU Time (total) 

32 

3.1 ■ 10-'= 

0.6369899mc 

5.8 

■ 10-2 

40.2 min 

5.3 ■ 102 hrs (est.) 



0.6369792z.« 

5.6 

■ 10-4 


6 sec 

64 

1.6- 10-2 

0.6367165mc 

1.5 

■ 10-4 

79.1 min 

1.0 ■ 104 hrs (est.) 



0.6367099l« 

1.4 

■ 10-4 


13 sec 

128 

7.8 ■ 10-2 

0.6366653mc 

7.2 

■ 10-2 

171 min 

2.2 ■ 104 hrs (est.) 



0.6366423l« 

3.5 

■ 10-2 


26 sec 

256 

3.9 ■ 10-2 

0.6366388mc 

3.0 

■ 10-2 

355 min 

4.7 ■ 104 hrs (est.) 



0.6366254i« 

8.9 

■ 10-2 


53 sec 

512 

2.0 ■ 10-2 

0.636621 8mc 

3.2 

■ 10-2 

705 min 

9.4 ■ 104 hrs (est.) 



0.6366212z.« 

2.2 

■ 10-2 


1.8 min 



0.6366198„ac/ 






certain cases EHia, our method is more universal one and is applicable to a wide class of potentials as a unified 
approach. It needs a constant memory size, which depends only on the final mesh size M and the rank r of the matrix 
of solution from ( |30l ) at each iteration step. Its complexity, then, is similar to the classical time-stepping schemes for 
the solution of the reaction-diffusion equations in a bounded domain. It also shows a favourable scaling. 

It is natural to extend the approach presented in the current work to higher dimensions. Then, instead of one¬ 
dimensional convolutions we will have to work with (f-dimensional convolutions, where d is the dimension of the 
problem. The extended domain will be [-na,na]‘^, where n is the number of time steps (equal to the dimension of 
the path integral). Thus, for higher dimensions the solution can be treated as a (t/ H- l)-dimensional tensor of size 
M X n X ... X n. Instead of the matrix low-rank approximation, stable low-rank factorization based on the tensor 
train decomposition Il44l could be used, with the final cost approximately equal to the cost of the computation the 
convolutions on the small domain. 

Finally, the most intriguing part of the work to be done is to apply the similar techniques to the Schrodinger 
equation. There, the convolution is no longer a convolution with a Gaussian function. Thus, the problem is much 
more difficult and our approach requires modifications. The presented method can also be applied to path integrals 
arising in other application areas, including the financial mathematics. The main requirement is that the integrand 
depends on the sum of variables multiplied by a separable function. 
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Appendix A. The cross approximation of matrices 

Let A 6 and I - {ii,i 2 , L), J - { 71,72 ■ • ■ 7 V). be subsets of / = {!,...,«) and J - {1,... ,m), respectively, 
r < min(n, m). By A = A(/, J) we denote a submatrix of A formed by the entries of A at the intersections of rows 
i e 1 and columns 7 e J. In this paper we use the following concept of the skeleton decomposition ll94l 1^ 1^ l97l . 
For any matrix A e of rank r there exist its decomposition 

A = BA^*C^, (A.l) 

where B = A(/, J), = A(/, J), and A = A(/, J) 6 is a certain submatrix of A, such that detA + 0. For the 

numerical reasons it is more effictive to work with orthogonal matrices. The decomposition (|A.l|i can be rewritten by 
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the factorization of the matrices B = QbRb and QR-decomposition, and by further factorization 

of the rank-r square matrix by the singular value decomposition (SVD) Il98ll99l . Thus, we will 

use the dyadic representation of ( |A.1[ ) 

A = XY^ X = QBUA2:y^ (^2) 

9=1 

If the rank of the matrix A is greater then r, in practice instead of exact equation ( |A.l| i we consider approximation 
in some norm. To obtain the decomposition ( |A.2[ ) in this case we use the cross approximation algorithm II1001IIOII 
based on the concept of the maximum volume submatrix (maxvol) introduced in II102111031 . We have implemented 
our version of the algorithm available at II104II1051 . Example of the usage of our code can be found in II106I . The 
new version of the code for complex and real matrices will be aviable soon at II107I . 

The rank in the cross approximation technique is determined adaptively. The algorithm starts from the guess 
rank r^ and at each iteration step k the subspace of vectors of B and is doubled (they are chosen by the maxvol 
subroutine, which returns a set of 2r^. row (column) indeces of a submatrix of (almost) maximum volume). The next 
value of the rank r^+i, r^+x < 2rk is chosen from the singular values of the matrix Xa of size 2 rt x 2 rk according to the 
following criterion 


rk+x = min {s | ^(s) < ed , 

l<5<2ri- 




A 


y2r/. 2 

^i=s+\ ^ i 

t2x^ 


= 0 . 


(A.3) 


The algorithm stops when ||X® - < Ec^a^** for the relative accuracy Sc- 

Approximation ( |A.2| i can be obtained by the SVD decomposition of the whole matrix A with 0{{n^ + m^)m) com¬ 
plexity, which is prohibitively slow. In contrast, the rank-r cross approximation requires only 0{{n + m)r) evaluations 
of the elements and 0((n + m)r^) additional operations. This becomes cmcial in practice, when the matrix element 
Aij is a time-consuming function to be calculated in a point (/, j) for a finite time or the given matrix is very large. 


Existence of such an approximation and convergence of the cross algorithm are discussed in Section 3.4 and Appendix 


Appendix B. Numerical investigation of the low-rank structure of the solution basis set 

Suppose, a function t{x) e Li)®) c™ be expanded into a series 

OO 

Kx) = ^ C/0/(x), Cl = 

1=0 


f 

kJ—O 


t(x) ipi(x)dx. 


(B.l) 


where 


4>i(x) 


1 






Hq{x) = 1 , H[{x) = 2x, H2(x) = 4x^ — 2 , 


\in\sjhl 

and Hi{x), are Hermite polynomials, with fast decaying coefficients ci, such that for a given accuracy ei 


(B.2) 


3^0 -xik) < ei;r(0)> 


X(l') = Xi 

1 = 1 '+ 1 


And the approximated function 


^0 

?(x) = ^ Q0/(x), ||f(x) - f(x)|| < Si 

/=! 


is of a canonical ei-rank Iq. Then, the question is about the rank structure of the matrix O/ constmcted as a reshape of 
a corresponding one-dimensional basis vector 0/(x,) defined on the uniform mesh (ffS]) 
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(B.3) 


















If the matrices are the low-rank ones then the matrix of a target function t(x) 


T = ^ c/ O/, Tij^ t(xi+j.M), 


(BA) 


1=1 


is of low-rank as well, which does not exceeds the upper bound Iq ■ rmax, where r^ax = maxi<;<;„(rank(0;)), but 


practically, it is of order r^ax- In the Table B.4 we present first several singular values of matrix O;. As it can be seen, 
each matrix has the low-rank structure. It would be nice to prove this numerical fact theoretically. 

Finally, to estimate the rank of the approximation ( |B.4| i one needs only to compute the coefficients c/ in the 
expansion ( |B.1[ ) and investigate their behaviour. This idea is similar to the QTT approach applied to the Laplace and 
its inverse operators in lllOSi . 


Appendix C. Proof of the lemmas 

Proof. (Lemmaj^. For the basis set the following equality holds 


n-i 


I/tt “ 'y ] 


i=0 


Then, according to (|^ 


/=0 


/=0 


Proof. (Lemma|^. From the equality r,„ = i 


it follows that 




r2-l 


/=0 


r2-l 


r2-l 


(c.i) 


ri~l ri-I ri-1 

^ ttm/uf , 0 = ^ ami [uf, O]^ , O \.m^^ £r,„;U;, Vm G [0,k - 1]. (C.2) 

J/r i'=o 


□ 


0^, , « R,„ = Vm G [0, k - 1]. 

H i=0 


i=0 


□ 


Proof. (Lemmaj^. From definition ( |22| i and the decomposition Im = Z-lo' Tm/U/, Sjm G [0,k], i 


it follows that 




ri-1 

“ Im “ ^ ^ ymi^i ~ ^ ^ Tmi 
/-O /-o 


n-1 


W; 

Um,(M-\) 


n-1 


fff7 “ 'y ] ymi^i- 
i=0 


□ 


Appendix D. Solution for the harmonic oscillator potential 

In this section we analytically integrate equations ([^, (|^, ( [T0| l with the initial condition ( |3T| l. 
Let us define F^"'(x) s 'B^^\x) for harmonic potential (|3T|l. Starting from k - n. 


'Ff (x) = 





and making use of the integral 

X oo _ 

p(J3,y + f)p(a,f)df^p[-^,y)^ y G M, 

oo ~ 

16 


a,/3 > 0, 













Table B.4: Singular values (s.v.) of the matrix ^B.3^ composed from the discretized basis jB.2) . The order of polynomials is labeled by /, (T\ is the 
first (absolute) singular value, then (Tij(T\ corresponds to the following relative singular values. The size of the matrix is 8000 x 1024. As it can be 
seen (numerical) ranks does not exceed the value of 8 and grow from small to bigger /. 


1 

0-1 

0 - 2 / 0-1 

0 - 3 / 0 -! 

0 - 4 / 0 -! 

0 - 5 / o-i 


0 - 6 / o-i 


0 - 7 / 0 -! 

o ' s / o'i 

0 - 9 / o-i 

0 

32.2 

0.96 

9.0 ■ 10-^ 

8.7 ■ 10-3 

8.9- 

10- 

16 

8.7- 

10- 

16 

7.4 

10- 

16 

6.4 

10- 

16 

5.5 

10- 

16 

1 

35.5 

0.77 

0.6 ■ 10-3 

0.6 ■ 10-3 

1.1 ■ 

10- 

14 

1.1 ■ 

10- 

14 

3.7 

10- 

-16 

3.4 

10- 

-16 

2.8 

10- 

-16 

2 

40.3 

0.48 

2.2 ■ 10-3 

1.6-10-3 

8.6- 

10- 

14 

8.3 ■ 

10- 

14 

5.6 

10- 

-16 

5.0 

10- 

-16 

4.6 

10- 

-16 

3 

38.0 

0.62 

0.8 ■ 10-3 

0.7 ■ 10-3 

6.2- 

10- 

13 

5.1 ■ 

10- 

13 

5.0 

10- 

-16 

4.4 

10- 

-16 

4.4 

10- 

-16 

4 

37.0 

0.68 

2.1 ■ 10-3 

1.8 ■ 10-3 

3.9- 

10- 

12 

3.8 ■ 

10- 

12 

6.6 

10- 

-16 

4.6 

10- 

-16 

4.5 

10- 

-16 

5 

35.5 

0.77 

5.0 ■ 10-3 

4.1 ■ 10-3 

2.1 ■ 

10- 

11 

1.5 ■ 

10- 

11 

5.0 

10- 

-16 

4.0 

10- 

-16 

4.0 

10- 

-16 

6 

38.3 

0.59 

9.2 ■ 10-3 

5.8 ■ 10-3 

9.2- 

10- 

11 

8.8 ■ 

10- 

11 

5.3 

10- 

-16 

4.8 

10- 

-16 

4.6 

10- 

-16 

7 

33.9 

0.83 

0.18 

0.13 

4.8 ■ 

10- 

10 

3.7 ■ 

10- 

10 

3.8 

10- 

-16 

3.8 

10- 

-16 

3.4 

10- 

-16 

8 

36.2 

0.68 

0.25 

0.09 

1.9 

10- 

-9 

1.8 

10- 

-9 

4.1 

10- 

-16 

4.0 

10- 

-16 

3.5 

10- 

-16 

9 

38.2 

0.47 

0.29 

0.26 

7.3 

10 

-9 

7.0 

10- 

-9 

4.7 

10- 

-16 

3.8 

10- 

-16 

3.6 

10- 

-16 

10 

28.7 

0.85 

0.60 

0.58 

3.6 

10' 

-8 

2.9 

10- 

-8 

7.9 

10- 

-16 

4.7 

10- 

-16 

4.5 

10- 

-16 

11 

31.1 

0.73 

0.65 

0.35 

1.1 

10- 

-7 

9.4 

10- 

-8 

5.5 

10- 

-16 

5.0 

10- 

-16 

4.7 

10- 

-16 

12 

33.1 

0.65 

0.58 

0.25 

3.7 

10- 

-7 

3.1 

10- 

-7 

5.3 

10- 

-16 

5.0 

10- 

-16 

4.8 

10- 

-16 

13 

30.1 

0.69 

0.68 

0.51 

1.3 

10 

-6 

1.2 

10 

-6 

7.0 

10- 

-16 

5.7 

10- 

-16 

5.3 

10- 

-16 

14 

26.2 

0.95 

0.74 

0.67 

3.9 

10 

-6 

3.2 

10- 

-6 

7.2 

10- 

-16 

6.9 

10- 

-16 

6.1 

10- 

-16 

15 

30.2 

0.74 

0.71 

0.38 

8.0 

10 

-6 

8.0 

10- 

-6 

6.2 

10- 

-16 

6.0 

10- 

-16 

5.3 

10- 

-16 

16 

30.8 

0.82 

0.63 

0.21 

2.0 

10- 

-5 

1.7 

10- 

-5 

5.3 

10- 

-16 

4.9 

10- 

-16 

4.6 

10- 

-16 

17 

28.4 

0.92 

0.67 

0.43 

5.5 

10- 

-5 

5.2 

10- 

-5 

5.1 

10- 

-16 

4.6 

10- 

-16 

4.6 

10- 

-16 

18 

28.5 

0.88 

0.68 

0.48 

1.3 

10- 

-4 

1.1 

10- 

-4 

5.9 

10- 

-16 

5.5 

10- 

-16 

5.4 

10- 

-16 

19 

30.0 

0.86 

0.51 

0.47 

2.6 

10- 

-4 

2.3 

10 

-4 

5.0 

10- 

-16 

4.8 

10- 

-16 

4.5 

10- 

-16 

20 

31.6 

0.80 

0.56 

0.24 

5.3 

10- 

-4 

3.4 

10- 

-4 

5.3 

10- 

-16 

4.7 

10- 

-16 

4.5 

10- 

-16 

21 

28.9 

0.90 

0.72 

0.23 

1.12 

■ 10 

-3 

1.1 

10- 

-3 

4.6 

10- 

-16 

4.4 

10- 

-16 

4.3 

10- 

-16 

22 

27.0 

0.98 

0.76 

0.47 

2.5 

10 

-3 

2.4 

10 

-3 

6.2 

10- 

-16 

5.6 

10- 

-16 

4.8 

10- 

-16 

23 

30.3 

0.67 

0.62 

0.59 

4.03 
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where p(a, x) is defined in ([T7]i, we obtain 
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By induction, we conclude that 


\\iin)f \ r^(n) I^k-l n ‘^yk ^ 

—e, Pk-i^- -, yk^Pk + wk 

V TT d + yj: 


6 t 


1 + (n - k)5t 


p(n) _ IA Pn-l 

^ V Tn V T«-I 


/—, \ <k <n 


References 

1. Feynman, R.P.. Space-time approach to non-relativistic quantum mechanics. Rev Mod Phys 1948;20:367-387. URL: http: //link. aps. 
org/doi/10.1103/RevModPhys.20.367 doi 10.1103/RevModPhys.20.367 

2. Feynman, R., Hibbs, A.. Quantum mechanics and path integrals. International series in pure and applied physics; McGraw-Hill; 1965. 

3. Garrod, C.. Hamiltonian path-integral methods. Rev Mod Phys 1966;38:483^94. doi 10.1103/RevModPhys. 38.483 

4. Abrikosov, A.A., Dzyaloshinskii, I., Gorkov, L.R. Methods of quantum field theory in statistical physics. Dover; 1975. URL: https: 
//cds.cern.ch/record/107441 

5. Mahan, G.. Many-Particle Physics. Physics of Solids and Liquids; Springer; 2000. URL: http://books, google. co.uk/books?id= 
xzSgZ4-yyMEC 

6. Norman Bleistein, R.A.H.. Asymptotic Expansions of Integrals. Dover Edition; 2010. 

7. Zinn-Justin, J.. Quantum field theory and critical phenomena (3d edition). Clarendon Press, Oxford; 1996. 

8. Polonyi, J.. Lectures on the functional renormalization group method. Central European Journal of Physics 2003',\i].). URL: http: 
//dx.doi.org/10.2478/BF02475552 doi 10.2478/BF02475552 

9. Dick, J., Kuo, F., Peters, G., Sloan, I.. Monte Carlo and Quasi-Monte Carlo Methods. Springer; 2012. doi 10.1007/ 
978-3-642-41095-6 

10. Holtz, M.. Sparse Grid Quadrature in High Dimensions with Applications in Finance and Insurance. 2011. doi 10.1007/ 
978-3-642-16004-2 

11. Garcke, J., Griebel, M.. Sparse Grids and Applications. Springer; 2013. doi 10.1007/978-3-642-31703-3 

12. Wong, K.Y.. Review of feynmans path integral in quantum statistics: from the molecular schrodinger equation to kleinerts variational 
perturbation theory. Communications in Computational Physics 2014; 15(4). doi 10.4208/cicp. 140313.070513s 

13. Masujima, M.. Path Integral Quantization and Stochastic Quantization. Springer; 2008. 

14. Kleinert, H.. Path Integrals in Quantum Mechanics, Statistics, Polymer Physics, and Financial Markets, 5th ed. World Scientific; 2009. 

15. Crank, J.. The Mathematics of Diffusion. Clarendon press, Oxford; 1975. 

16. Bass, R.F.. Diffusions and Elliptic Operators. Springer; 1998. 

17. Chaichian, A., Demichev, A.. Path Integrals in Physics; vol. I. loP publishing; 2001. 

18. Borodin, A.N., Salminen, P. Handbook of Brownian motion. Facts and formulae. Springer; 2002. 

19. Kac, M.. On some connections between probability theory and differential and integral equations. 1951. URL: http://projecteuclid. 
org/euclid.bsmsp/1200500229 

20. Karatzas, L, Shreve, S.E.. Brownian motion and partial differential equations. Springer; 1991. 

21. Mazo, R.M.. Brownian motion. Fluctuations, Dynamics, and Applications. Clarendon Press, Oxford; 2002. 

22. Nelson, E.. Dynamical Theories of Brownian motion. Princeton University Press; 2001. 

23. Smolyak, S.A.. Quadrature and interpolation formulas for tensor products of certain class of functions. Dokl Akad Nauk SSSR 
1964;148(5):1042-1053. Transl.: Soviet Math. Dokl. 4:240-243, 1963. 

24. Gerstner, T., Griebel, M.. Dimension-adaptive tensor-product quadrature. Computing 2003;71:65-87. doi 10.1007/ 
S00607-003-0015-5 

25. Gerstner, T., Griebel, M.. Numerical integration using sparse grids. Numerical Algorithms l99SA^{3-4):209-232. URL: http://dx. 
doi.org/10.1023/A%3A1019129717644 doi: 10.1023/A:1019129717644 


18 













26. Makri, N., Miller, W.H.. Monte carlo integration with oscillatory integrands: implications for feynman path integration in real time. Chem¬ 
ical Physics Letters 1987;139(1):10 - 14. URL: http://www.sciencedirect.coin/science/article/pii/0009261487801422 
doi http://dx.doi.org/10.1016/0009-2614(87)80142-2 

27. Gorshkov, V.N., Tretiak, S., Mozyrsky, D.. Semiclassical monte-carlo approach for modelling non-adiabatic dynamics in extended 
molecules. Nat Commun 2013;4:1 8. doi 10.1038/nconims3144 

28. de Lathauwer, L.. A survey of tensor methods. In: IEEE International Symposium on Circuits and Systems. 2009:2113-2116. doi 10. 
1109/iscas.2009.5118377 

29. Smilde, A., Bro, R., Geladi, R. Multi-way analysis with applications in the chemical sciences. Wiley; 2004. 

30. Khoromskij, B.N.. Introduction to tensor numerical methods in scientific computing. Preprint, Lecture Notes 06-2011; University of 
Zurich; 2010. URL: http://www.math.uzh. ch/fileadmin/math/preprints/06_ll .pdf 

31. Khoromskij, B.N.. Tensor-stmctured numerical methods in scientific computing: Survey on recent advances. Chemometr Intel! Lab Syst 
2012;110(1):1-19. doi 10.1016/j .chemolab.2011.09.001 

32. Hitchcock, F.L.. The expression of a tensor or a polyadic as a sum of products. J Math Phys 1927;6(1):164-189. 

33. Hitchcock, F.L.. Multiple invariants and generalized rank of a p-way matrix or tensor. J Math Phys 1927;7(l):39-79. 

34. Grasedyck, L., Kressner, D., Tobler, C.. A literature survey of low-rank tensor approximation techniques. GAMM-Mitteilungen 
2013;36(l):53-78. URL: http://dx.doi.org/10.1002/gainm.201310004 doi 10.1002/gainm.201310004 

35. Kolda, T.G., Bader, B.W.. Tensor decompositions and applications. SIAM Review 2009',5l{3):455-500. doi 10.1137/07070111X 

36. Grasedyck, L., Wolfgang, H.. An introduction to hieraixhical h-rank and tt-rank of tensors with examples. Computational Methods in 
Applied Mathematics 201 l;n(3):291. 

37. Oseledets, I.V., Savostianov, D.V., Tyrtyshnikov, E.E.. Tucker dimensionality reduction of three-dimensional arrays in linear time. S/AM 
J Matrix Anal Appl 200H;30{3) :939-956. doi 10.1137/060655894 

38. Tucker, L.R.. Some mathematical notes on three-mode factor analysis. Psychornetrika 1966;31:279-311. doi 10.1007/BF02289464 

39. Tucker, L.. Implications of factor analysis of three-way matrices for measurement of change. Problems in measuring change 1963;: 122- 
137. 

40. Tucker, L.R.. The extension of factor analysis to three-dimensional matrices. Contributions to mathematical psychology 1964;: 109-127. 

41. de Lathauwer, L.,deMoor, B., Vandewalle, J.. Amultilinearsingularvaluedecomposition. S/AM/MflrrijcAn(3/App/2000;21:1253-1278. 
doi 10.1137/S0895479896305696 

42. Oseledets, I.V., Tyrtyshnikov, E.E.. TT-cross approximation for multidimensional arrays. Linear Algebra Appl 20\0'A32{1):10-SS. 
doi 10.1016/j.laa.2009.07.024 

43. Oseledets, I.V., Tyrtyshnikov, E.E.. Breaking the curse of dimensionality, or how to use SVD in many dimensions. SIAM J Sci Comput 
2009;31(5):374^3759. doi 10.1137/090748330 

44. Oseledets, I.V.. Tensor-train decomposition. SIAM J Sci Comput 20\\\33{5):2295-23\l. doi 10.1137/090752286 

45. Hackbusch, W., Kiihn, S.. A new scheme for the tensor representation. J Fourier Anal Appl 2009',\5{5):106-122. doi 10.1007/ 
S00041-009-9094-9 

46. Grasedyck, L.. Hierarchical singular value decomposition of tensors. SIAM J Matrix Anal Appl 20\0',3\{A):2.029—2054. doi 10.1137/ 
090764189 

47. Hackbusch, W.. Tensor spaces and numerical tensor calculus. Springer-Verlag, Berlin; 2012. ISBN 978-3642280269. 

48. Khoromskij, B.N., Khoromskaia, V., Chinnamsetty, S.R., Flad, H.J.. Tensor decomposition in electronic structure calculations on 3D 
Cartesian grids. J Comput Phys 2009',22%{\6):5149-5162. doi 10.1016/j . jcp. 2009.04.043 

49. Khoromskij, B.N., Khoromskaia, V.. Multigrid accelerated tensor approximation of function related multidimensional aiTays. SIAM J Sci 
Comput 2009;31(4):3002-3026. doi 10.1137/080730408 

50. Khoromskaia, V., Khoromskij, B.. Tensor numerical methods in quantum chemistry: from hartree-fock to excitation energies. Phys Chem 
Chem Phys 2015;URL: http: //arxiv. org/abs/1504.06289 doi: 10.1039/c5cp01215e 

51. Schneider, R., Rohwedder, T., Blauert, J.. Direct minimization for calculating invariant subspaces in density functional computations of 
the electronic structure. Journal of Comp Math 2009;27:360. 

52. Flad, Heinz-Jrgen,, Hackbusch, Wolfgang,, Schneider, Reinhold,. Best n-term approximation in electronic structure calculations, ii.jastrow 
factors. ESAIM: M2AA 2007;41(2):261-279. URL: http://dx.doi.Org/10.1051/in2aii:2007016 doi 10.1051/in2an:2007016 

53. Khoromskij, B.N., Khoromskaia, V., Flad., H.J.. Numerical solution of the Hartree-Fock equation in multilevel tensor-structured format. 
SIAM J Sci Comput 201 V,33{l):45-65. doi 10.1137/090777372 

54. Kazeev, V.A., Khoromskij, B.N.. Low-Rank explicit QTT representation of the Laplace operator and inverse. SIAM journal on matrix anal¬ 
ysis and applications 2012;33(3):742-758. URL: http://www.inis.mpg.de/de/publicatioiis/preprints/2010/prepr2010-75. 
html doi 10.1137/100820479 

55. Hackbusch, W., Khoromskij, B.N.. Low-rank Kronecker-product approximation to multi-dimensional nonlocal operators. I. Separable 
approximation of multi-variate functions. Computing 2006;76(3-4):177-202. doi: 10.1007/s00607-005-0144“0 

56. Hackbusch, W., Khoromskij, B.N.. Low-rank Kronecker-product approximation to multi-dimensional nonlocal operators. II. HKT repre¬ 
sentation of certain operators. Computing 2006',16{3-4):203-225. doi 10.1007/s00607-005-0145-z 

57. Beylkin, G., Mohlenkamp, M.J.. Numerical operator calculus in higher dimensions. Proc Nat Acad Sci USA 2002;99{16):10246-\0251. 
doi 10.1073/pnas.112329799 

58. Gavrilyuk, I.R, Hackbusch, W., Khoromskij, B.N.. Tensor-product approximation to the inverse and related operators in high-dimensional 
elliptic problems. Computing 2005;(74):131—157. 

59. Khoromskij, B.N.. Tensor-structured preconditioners and approximate inverse of elliptic operators in M.^. Constr Approx 2009;(30):599- 
620. doi 10.1007/S00365-009-9068-9 

60. Khoromskij, B.N., Khoromskaia, V.. Low rank Tucker-type tensor approximation to classical potentials. Central European journal of 
mathematics 2001 \5{3):523-550. doi 10.2478/sll533-007-0018-0 

61. Hackbusch, W., Khoromskij, B.N.. Tensor-product approximation to operators and functions in high dimensions. Journal of Complexity 


19 


2007;23(46):697 - 714. URL: http://www.sciencedirect.coin/science/article/pii/S0885064X07000532 doi http://dx. 
doi. org/10.1016/j . jco. 2007.03.007, festschrift for the 60th Birthday of Henryk WoAniakowski. 

62. Hackbusch, W., Braess, D.. Approximation of ^ by exponential sums in [1, oo]. IMA JNumer Anal 2005;25(4):685-697. 

63. Oseledets, I.V., Muravieva, E.A.. Fast orthogonalization to the kernel of discrete gradient operator with application to the Stokes problem. 
Linear Algebra Appl 2010;432(6): 1492-1500. doi 10.1016/j . laa. 2009.11.010 

64. Oseledets, I.V.. Constructive representation of functions in low-rank tensor formats. Constr Appr URL: http://pub. 

inni.ras.ru/pub/ininras2010-04.pdf doi 10.1007/s00365-012-9175-x 

65. Beylkin, G., Mohlenkamp, M.J.. Algorithms for numerical analysis in high dimensions. SIAM J Sci Comput 2005;26{6)\2133-2159. 

66. Dolgov, S.V., Khoromskij, B.N., Oseledets, I.V.. Fast solution of multi-dimensional parabolic problems in the tensor train/quantized 
tensor train-format with initial application to the Fokker-Planck equation. SIAM J Sci Comput 2012;34(6):A3016-A3038. doi 10.1137/ 
120864210 

67. Khoromskij, B.N., Oseledets, I.V.. Quantics-TT collocation approximation of parameter-dependent and stochastic elliptic PDFs. Comput 
Meth Appl Math 20l0;W{4):316-394. doi 10.2478/cinani-2010-0023 

68. Khoromskij, B.N., Oseledets, I.V.. DMRG+QTT approach to computation of the ground state for the molecular Schrodinger operator. 
Preprint 69; MPI MIS; Leipzig; 2010. URL: http://www.mis .mpg.de/preprints/2010/preprint2010_69.pdf 

69. Oseledets, I.V.. Approximation of 2^ x 2^ matrices using tensor decomposition. SIAM J Matrix Anal Appl 2010;31(4):2130-2145. 
doi 10.1137/090757861 

70. Lebedeva, O.S.. Tensor conjugate-gradient-type method for Rayleigh quotient minimization in block QTT-format. Russ J Numer Anal 
Math Modelling 2Qi\\\26{5yA654'&9. doi 10.1515/r jnamm. 2011.026 

71. Khoromskij, B.N., Oseledets, I.V.. QTT-approximation of elliptic solution operators in high dimensions. Rus J Numer Anal Math Model 
2011;26(3):303-322. doi 10.1515/r jnamm. 2011.017 

72. Savostyanov, D.V.. QTT-rank-one vectors with QTT-rank-one and full-rank Fourier images. Preprint 45; MPI MIS; Leipzig; 2011. URL: 
http://www.mis.mpg.de/preprints/201l/preprint2011_45.pdf 

73. Khoromskij, B.N.. logn)-Quantics approximation of N-d tensors in high-dimensional numerical modeling. Constr Appr 

2011;34(2):257-280. doi 10.1007/s00365-011-9131-l 

74. Khoromskaia, V., Khoromskij, B.N., Schneider, R.. Tensor-structured factorized calculation of two-electron integrals in a general 
basis. SIAM Journal on Scientific Computing 2013;35(2):A987-A1010. URL: http://dx.doi.org/10.1137/120884067 doi 10. 
1137/120884067 arXiv:http://dx.doi.org/10.1137/120884067 

75. Ballani, J.. Fast evaluation of singular bem integrals based on tensor approximations. Numerische Mathematik 20i2;121(3):433—460. URL: 
http://dx.doi.Org/10.1007/s00211-011-0436-6 doi 10.1007/s00211-011-0436-6 

76. Rakhuba, M.V., Oseledets, I.V.. Fast multidimensional convolution in low-rank tensor formats via cross approximation. SIAM Jour¬ 
nal on Scientific Computing 2015;37(2):A565-A582. URL: http://dx.doi.org/10.1137/140958529 doi 10.1137/140958529 
arXiv:http://dx.doi.org/10.1137/140958529 

77. Khoromskij, B.N.. Fast and accurate tensor approximation of multivariate convolution with linear scaling in dimension. J Comp Appl Math 
2010;234(11):3122-3139. doi: 10.1016/j . cam. 2010.02.004 

78. Kazeev, V., Khoromskij, B.N., Tyrtyshnikov, E.E.. Multilevel Toeplitz matrices generated by tensor-structured vectors and convolution 
with logarithmic complexity. Tech. Rep. 36; MPI MIS; Leipzig; 2011. URL: http://www.mis.mpg.de/publications/preprints/ 
201l/prepr2011-36.html 

79. Dolgov, S.V., Khoromskij, B.N., Savostyanov, D.V.. Superfast Fourier transform using QTT approximation. J Fourier Anal Appl 
2012;18(5):915-953. doi 10.1007/s00041-012-9227-4 

80. Beylkin, G., Monzon, L.. Approximation by exponential sums revisited. Appl Comput Harm Anal 2010;28(2):131-149. doi 10.1016/j . 
acha.2009.08.011 

81. Brigham, E.O.. The fast Fourier transform and its applications. Prentice Hall; 1988. 

82. Nussbaumer, H.J.. Fast Fourier transform and convolution algorithms. Springer; 1981. 

83. Chiu, J., Demanet, L.. Sublinear randomized algorithms for skeleton decompositions. SIAM Journal on Matrix 
Analysis and Applications 2013;34(3):1361-1383. URL: http://dx.doi.org/10.1137/110852310 doi 10.1137/110852310 
arXiv:http://dx.doi.org/10.1137/110852310 

84. Bebendorf, M., Grzhibovskis, R.. Accelerating galerkin bem for linear elasticity using adaptive cross approximation. Mathematical 
Methods in the Applied Sciences 2006;29(14): 1721-1747. URL: http: //dx. doi. org/10.1002/mma. 759 doi: 10.1002/mma. 759 

85. Bebendorf, M., Kunis, S.. Recompression techniques for adaptive cross approximation. J Integral Equations Applications 2009;21{3):33l- 
357. URL: http://dx.doi.org/10.1216/JIE-2009-21-3-331 doi 10.1216/JIE-2009-21-3-331 

86. Bebendorf, M.. Adaptive cross approximation of multivariate functions. Constructive Approximation 201 ].;34(2)\\49-ll9. URL: http: 
//dx.doi.org/10.1007/s00365-010-9103-x doi 10.1007/s00365-010-9103-x 

87. Hairer, E., Lubich, C., Wanner, G.. Geometric Numerical Integration. Structure-Preserving Algorithms for Ordinary Differential Equations. 
Springer series in computational mathematics; Springer; 2006. doi 10.1007/3-540-30666-8 

88. Brezinski, C., Zaglia, M.R.. Extrapolation methods. Theory and practice. Elsevier; 1991. 

89. Stoer, J., Bulirsch, R.. Introduction to Numerical Analysis. 3rd ed. ed.; Springer; 2002. 

90. Makri, N.. Numerical path integral techniques for long time dynamics of quantum dissipative systems. Journal of Mathematical Physics 
1995;36(5):2430-2457. doi http: //dx. doi. org/10.1063/1.531046 

91. Makri, N.. Blip decomposition of the path integral: Exponential acceleration of real-time calculations on quantum dissipative systems. The 
Journal of Chemical Physics 2014; 141 (13): 134117. doi http: //dx. doi. org/10.1063/1.4896736 

92. Dubach, E.. Artificial boundary conditions for diffusion equations: Numerical study. Journal of Computational and Applied Mathematics 
1996;70(1):127 - 144. URL: http: //www. sciencedirect. coin/science/article/pii/0377042795001409 doi http: //dx.doi. 
org/10.1016/0377-0427(95)00140-9 

93. Wu, X., Sun, Z.Z.. Convergence of difference scheme for heat equation in unbounded domains using aitificial boundary condi- 


20 


tions. Applied Numerical Mathematics 2004;50(2):261 - 277. URL: http://www.scieiicedirect.com/science/article/pii/ 
S0168927404000030 doi http://dx.doi.org/10.1016/j.apnnra.2004.01.001 

94. Tyrtyshnikov, E.E.. Mosaic-skeleton approximations. Calcolo \99&,'i'i{\yAl-51. doi: 10.1007/BF02575706 

95. Goreinov, S.A., Tyrtyshnikov, E.E., Zamarashkin, N.L.. Atheory of pseudo-skeleton approximations. I997;261:l- 

21. doi 10.1016/30024-3795(96)00301-1 

96. Goreinov, S.A., Tyrtyshnikov, E.E., Zamarashkin, N.L.. Pseudo-skeleton approximations of matrices. Reports of Russian Academy of 
Sciences 1995;342(2):151-152. 

97. Goreinov, S.A., Zamarashkin, N.L., Tyrtyshnikov, E.E.. Pseudo-skeleton approximations by matrices of maximum volume. 

Notes 1997;62(4):515-519. doi 10.1007/BF02358985 

98. Golub, G.H., Loan, C.F.V.. Matrix computations, 4th edition. 2012. 

99. Demmel, J.W.. Applied Numerical Linear Algebra. 1997. 

100 . Tyrtyshnikov, E.E.. Incomplete cross approximation in the mosaic-skeleton method. Computing 200Q;64{4y361-3S0. doi 10 . 1007 / 
S006070070031 

101. Bebendorf, M.. Approximation of boundary element matrices. Numer Mathem200036{4y.565-5^9. doi 10.1007/pl00005410 

102. Goreinov, S.A., Tyrtyshnikov, E.E.. The maximal-volume concept in approximation by low-rank matrices. Contemporary Mathematics 
2001;280:47-51. 

103. Goreinov, S.A., Oseledets, I.V., Savostyanov, D.V., Tyrtyshnikov, E.E., Zamarashkin, N.L.. How to find a good submatrix. Research 
Report 08-10; ICM HKBU; Kowloon Tong, Hong Kong; 2008. URL: http: //www.math.hkbu.edu.hk/ICM/pdf/08-10.pdf 

104. Litsarev, M.S., Oseledets, I.V.. Cross2d C++code: included in the DEPOSIT destribution. 2014. URL: https://bitbucket. org/ 
appl_m729/code-deposit 

105. Litsarev, M.S., Oseledets, I.V.. The DEPOSIT computer code based on the low rank approximations. Computer Physics Communications 
2014;185(10):28012802. URL: http://www.sciencedirect.com/science/article/pii/S0010465514002173 doi http://dx. 
doi.org/10.1016/j.cpc.2014.06.012 

106. Litsarev, M.S., Oseledets, I.V.. Fast low-rank approximations of multidimensional integrals in ion-atomic collisions modelling. Numerical 
Linear Algebra with Applications 2015;URL: http://dx.doi. org/10.1002/nla. 2008 doi 10.1002/nla. 2008 

107. Litsarev, M.S.. Cross2d C++code for the real and complex matrices. Template version. 2015. URL: https://bitbucket. org/appl_ 
m729/dzcross2d 

108. Gavrilyuk, I.P, Khoromskij, B.N.. Quantized-TT-Cayley transform to compute dynamics and spectrum of high-dimensional Hamil¬ 
tonians. Computational methods in applied mathematics 2011;ll(3):273-290. URL: http://www.mis.mpg.de/de/publications/ 
preprints/2011/prepr2011-43.html doi 10.2478/cmam-2011-0015 


21 


